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Abstract 

A method based on sector decomposition has been developed to calculate the double real 
radiation part of the process e + e~ — > 3jets at next-to-next-to-leading order. It is shown in an 
example that the numerical cancellation of soft and collinear poles works well. The method 
is flexible to include an arbitrary measurement function in the final Monte Carlo program, 
such that it allows to obtain differential distributions for different kinds of observables. This 
is demonstrated by showing 3-, 4- and 5-jet rates at order for a subpart of the process. 



1 Introduction 



Experiments at LEP have shown that the measurement of jet rates and shape observables in 
e + e~ collisions allow for very stringent tests of the Standard Model, in particular of predictions 
relying largely on Quantum Chromodynamics (QCD) [1], allowing for example a very precise 
determination of the strong coupling constant a s . A precise knowledge of a s in turn is of major 
importance at hadron colliders, especially at the LHC. However, the measurements from jets 
and shapes in e + e~ collisions, although being very precise, have not been included in the world 
average value for a s , because it is based only on measurements where next-to-next-to leading 
order (NNLO) theory predictions are available [2], while for e + e~ —>■ 3 jets, full NNLO predictions 
do not exist yet. 

A future International Linear Collider will allow for precision measurements at the per-millc 
level, which offer the possibility of a determination of a s with unprecedented precision. However, 
this will only be possible if the theoretical error can keep up with such a precision. As the 
present error on the NLO prediction for e + e~ — ► 3 jets is dominated by scale uncertainties [3-6], 
the calculation of the NNLO corrections to this process will surely lead to an important gain in 
precision. 

After the virtual two-loop corrections entering this calculation have become available [7-9], the 
bottleneck now is given by the real radiation part where up to two partons can become theoretically 
unresolved (soft and/or collinear), leading to infrared singularities upon phase space integration. 
These singularities have to be subtracted and cancelled with the ones from the virtual contributions 
before a Monte Carlo program can be constructed. At NNLO, the infrared singularities can be 
entangled in a complicated way, which renders the extraction of the poles a formidable task. Two 
different approaches can be followed to achieve this task: 

1 . Construction of a subtraction scheme where the subtraction terms are integrated analytically 
in D = 4 — 2e dimensions over the unresolved phase space, thus extracting the poles in 1/e. 
The main advantages of this approach are the following: It allows maximal (i.e. analytical) 
control over the pole terms, and it insures a minimal number of subtraction terms, as 
the latter are constructed manually by considering all physical situations where a singular 
configuration is approached. The drawbacks of this method are given by the fact that 
constructing such a scheme is a highly non-trivial and tedious task, especially in view of 
the fact that it is different for each colour structure. Further, the analytic integration over 
subtraction terms may become impossible when applying the method to other processes 
where several mass scales are involved. 

2. Sector decomposition, where the poles are isolated by an automated routine and the pole 
coefficients are integrated numerically. The advantages of this approach reside in the fact 
that the extraction of the infrared poles is algorithmic, being the same for all colour factors, 
and that the subtraction terms can be arbitrarily complicated as they are integrated only 
numerically. On the other hand, the algorithm which isolates the poles increases the number 
of original functions and in general does not lead to the minimal number of subtraction terms, 
thus producing rather large expressions. 

Approach 1. has been pursued by several groups in different variations [10-23], and the imple- 
mentation of the method based on antenna subtraction [19, 23] into a Monte Carlo program is 
presently underway [24]. The sector decomposition approach has seen a very rapid development 
recently. Sector decomposition is a general method to disentangle overlapping singularities in 



parameter space, originally used by K. Hepp [25] for overlapping ultraviolet singularities. It has 
been very successfully applied to various types of multi-loop integrals since [26-32]. Its applica- 
tion to NNLO phase space integrals, first proposed in [33], already lead to a number of promising 
results [15,32,34-37]. 

The present paper deals with the application of sector decomposition to the double real ra- 
diation part of the process e + e~ — > 3jets at NNLO, which involves subprocesses of the type 
7* — ► 5 partons. The matrix elements for the processes 7* — > qqggg and 7* — > qqq'q'g are huge, 
such that the calculation also involves non-trivial book-keeping and file handling tasks, which are 
not addressed in this article. The intention of this paper is to show that a method has been 
developed which can deal with 1 — > 5 parton processes efficiently, such that the construction of a 
fully differential Monte Carlo program for the process e + e~ — > 3 jets at NNLO is merely a matter 
of putting pieces together, which will be left to a future publication. Therefore, we only consider 
one sample topology (including the full tensor structure) as part of the full matrix element. For 
this topology, we first calculate the fully inclusive integral over the 5-parton phase space, leading 
to poles up to 1/e 4 . In order to prove the correctness of the result, we also calculate all possible 
cuts of this diagram with less than five particles in the final state. For the contribution from 
7* — > 4 partons, we use the result obtained in [35] by sector decomposition. The KLN theo- 
rem [38,39] guarantees that the sum of all possible cuts of the UV renormalised diagram is finite. 
This is demonstrated in section 2. However, the method is not limited to calculate only inclusive 
cross sections. As the singularities are disentangled by an algebraic algorithm, the inclusion of an 
arbitrary (infrared safe) measurement function - at the stage of the numerical evaluation of the 
finite functions produced by sector decomposition - does not present a problem. This is shown 
in section 3. As an illustration of the action of the measurement function, the 0(a^) 3~, 4- and 
5— jet rates are shown for the sample matrix element as a function of the cut parameter y within 
the JADE algorithm [40]. Section 4 contains the conclusions. Details of the calculation are given 
in appendix A. 

2 Cancellation of Divergences 

As explained above, the method presented here addresses the main difficulty in calculating the real 
radiation part of e + e~ — » 3jets at NNLO, which is the isolation and subtraction of the infrared 
poles which occur when integrating the squared amplitude over the phase space for 7* — > 5 partons. 

In order to check the correctness of the results for the integrals over the 1 — > 5 particle 
phase space, one can exploit the fact that the sum over all cuts of a given (UV renormalised) 
topology must be infrared finite. In order to demonstrate these cancellations, let us consider as 
an example the diagram depicted in Fig. Q occurring in the part ~ Cp of the squared amplitude 
for e+e" -> 3 jets at NNLO. 



Figure 1: The ladder diagram 




Summing over all cuts of this diagram and performing UV renormalisation, we obtain the 
condition 



Ti^ 5 + z\ Ti^ 4 + z 2 Ti^ 3 + z 3 Ti^ 2 = finite , 
where denotes the diagram with i cut lines as shown in Figure [21 
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Figure 2: Cancellation of IR divergences in the sum over all cuts of the renormalised graph 



2.1 UV renormalisation 

For i = 1,2, the renormalisation constants Zi (in Feynman gauge) already have been calculated 
in [35], to be given by 
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The 3-loop renormalisation constant Z3 will be derived in the following. Using the graphical 
BHPZ notation as in ref. [35], 3-loop renormalisation of the fermion selfenergy implies that the 
combination of graphs as shown in Fig. |3]is finite. 
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Figure 3: UV renormalisation of the quark propagator at C(af ) 
The explicit calculation yields: 
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The overall divergences = overall div of the diagrams above are thus given by 

+ 7 [| - y - I log (V/M 2 ) + \ log 2 (-p 2 /f< 2 )]} (2.8) 

+7 [ T - T - 5 1o s(-p 2 /" 2 ) + '°e 2 (-p 2 /m 2 )]} (2.9) 

3 



+\\\-^-\ lo S ("^/^) + \ log 2 ( V/V)] } (2.10) 



+- [3 - 4 " 3 log (V/V ) + I log 2 (V/V )]} • (2-11) 



e 6 

Note that we have adopted the MS prescription 
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From eqs. (j2.8j) to (|2.11|) we can now derive z% as 
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The non-local logarithmic terms cancel, as guaranteed by the BHPZ theorem [25,41,42]. 
2.2 Combining the renormalised diagrams 

In order to verify eq. ()2.1|) . we have to integrate the ladder diagrams corresponding to the process 
7* — > ipartons over the 1 — > i particle phase space. Up to i = 4, this has been done already in 
ref. [35], where Ti^4 has been calculated by sector decomposition. The important new ingredient 
here is the calculation of T%^. Before showing this calculation in more detail, let us first construct 



the expressions entering eq. (|2.1j) for i < 5. From ref. [35], we have 
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Combination of these results with the renormalisation constants given in section 12.11 leads to 

z 1 T 1 ^ 4 + z 2 T 1 ^ 3 + z 3 T 1 ^ 2 = (2.16) 
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What remains to be shown now is that the 5-parton contribution exactly cancels the poles 

in (t2~m 



2.3 Calculation of the 5-particle contribution 

The graph Ti_+5 is calculated numerically by sector decomposition. To this aim, the phase space 
integrals are brought to a form where all integrations are from zero to one, as described in more 
detail in appendix A. Note that the parametrisation given here is particularly convenient for the 
denominator structure of our sample topology. In order to deal with the full matrix element, 
several parametrisations have been worked out, each one optimised to be applied to a certain class 
of denominators. An automated subroutine scans the denominators of a given matrix element 
and applies the appropriate parametrisation. In this way, the full expression naturally is split into 
tractable subparts. 

After having performed the transformations of the phase space integration variables as ex- 
plained in appendix IA.11 the 1 — ► 5 phase space in D dimensions is given by ea. (|A.7[) : 

-1 10 

= 4 5 V) 2I? - 5 / ]\dt i [h(i-t s )]- i - e [t B (i-t 8 )t 10 (i-t 10 )]-^- e 

fa t 6 (l - te){l - t 7 )] 1_2e [(l - t 2 )t 3 {l - t 3 )t 4 (l - U)t 9 {l - t 9 )]- e t 2 7 ~ 3e (2.17) 



(5) 2^ 

r (4vr) 9 r(-2e)r(2 - 2e) ' 

Singularities only occur at the boundaries ti = 0,1. Further, one can split the integrations at 
ti = 1/2 and remap the variables to the unit cube to assure that all potential singularities occur 
only for ti — > 0. However, as this procedure doubles the number of integrals for each ti, it is 
only done for those variables where a singularity at t i = 1 is possible at all, in order to avoid a 
proliferation of terms. 

The matrix element typically contains terms of the structure 

I = [ dx [ dy x~ x ~ € (x + yy 1 , 
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where a naive subtraction of the singularity for x — > of the form 



dx / dyx~ ~ e f(x,y) = — / dy f(0,y)+ / dx dyx 
Jo e jo Jo Jo 



fails, because the singularities for x — ► and y — ► are overlapping. This is where sector 
decomposition shows its virtues. The working mechanism of sector decomposition already has 
been explained in detail in [27] and therefore will be outlined only shortly here. The basic idea is 
to first split the integration region into sectors where the variables x and y are ordered. 



J dx J dyx~ 1 ~ e (x + y)~ 1 [9(x-y) + @(y-x)] 



(1) (2) 

Then remapping the integration domain to the unit cube, the singularities in our simple example 
are already disentangled: After the substitutions y = xt in sector (1) and x = yt in sector (2), 
one has 

/ = f 1 dx x- l - e f 1 dt (1 + ty l + f 1 dy y- 1 -" ^ dt r 1-e (1 + t)- 1 . 
Jo Jo Jo Jo 

For more complicated functions, several iterations of this procedure may be necessary, but it is 
easily implemented into an automated subroutine. Once all singularities are factored out, the 
result can be expanded in e, where the subtraction of the pole terms naturally leads to plus 
distributions by the identity 
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In this way, a Laurent series in e is obtained, where the pole coefficients are sums of finite parameter 
integrals which can be evaluated numerically. 

Note that the numerator structure of the matrix element can only improve the infrared pole 
structure, such that it can be included later, at the stage of the expansion in e. It also should be 
mentioned that for some phase space parametrisations, required to tackle the full matrix element, 
square-root terms in the denominator are unavoidable. Such terms can spoil the simple scaling 
behaviour which is crucial for the algorithm to work. However, one can always find (nonlinear) 
variable transformations such that these terms can be mapped to a form which is amenable to 
sector decomposition. 

Applying the method to our sample diagram, and requiring a numerical precision of 1%, the 
following result is obtained after an integration time of about 20 minutes on a 2.8 GHz Pentium 
IV PC: 
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Combining eqs. (|2.16|) and (|2.18j) we see that all poles cancel within the numerical precision. 



3 Inclusion of a measurement function 

The isolation of infrared poles by sector decomposition is an algebraic procedure, leading to a set 
of finite functions for each pole coefficient as well as for the finite part. The finite part can be 
written to a Monte Carlo program and combined with any infrared safe measurement function. To 
this aim, one has to take the limit D — > 4 of the .D-dimensional phase space, which is non-trivial 
because in D = 4, the Gram determinant of five light-like momenta vanishes, which means that 
only 8 Mandelstam invariants are independent, whereas in D = 4 — 2e one has 9 independent 
invariants, i.e. 9 independent phase space integration variables, and sector decomposition acts 
in D = 4 — 2e dimensions. How this problem is solved is explained in appendix IA.21 It is also 
described there how the four-momenta of the particles in the final state in terms of energies and 
angles are reconstructed from the phase space integration variables ij. Note that the variables 
ti are transformed in the course of sector decomposition, such that for each function which is an 
endpoint of the sector decomposition tree, the expressions for the invariants Sy in terms of the final 
Monte Carlo integration variables look different. This requires careful (automated) book-keeping, 
but does not constitute a principal problem. 

Further, it has to be assured that the subtraction terms only come to action in phase space 
regions which are allowed by the measurement function. To illustrate this point, consider the 
simple one-dimensional example where the measurement function is just a step function Q(x — 
a), a > 0, and the "matrix element" after sector decomposition is given by a plus distribution 
[/(x)/x]+. If we naively combine the plus distribution with our measurement function, we obtain 



On the other hand, the /(0) term stems from the subtraction of a singularity at x = 0, which is 
now killed by our measurement function anyway, such that inclusion of the f(0) term would lead 
to a wrong result. Therefore, the correct way to include the measurement function is of course 



However, this does not mean that the e-expansions and subtractions have to be redone each time 
the measurement function is changed. It can be achieved by including symbolic functions in the 
e-expansion which, depending on the generated phase space point in the Monte Carlo program, 
take on the appropriate values. 

As an example, the JADE algorithm [40] to define 3-, 4- and 5— jet events has been implemented 
into a Monte Carlo program built upon the output of sector decomposition, using the multi- 
dimensional integration package BASES [43]. For the plots in Figs. and HJ), the diagram 
discussed in this paper (summed over all cuts) served as a toy matrix element, but it should 
be emphasised that the same Monte Carlo program can be used to calculate the full process 
e + e~ — ► 3, 4, 5 jets at order a^, once the contributions from the other topologies are implemented. 
Further, the architecture of the program is such that the JADE algorithm can be easily replaced 
by a different jet algorithm, and shape observables can also be defined. 

Figures and0J) show the 3-, 4- and 5-jet rates at order function of the jet separation 

parameter y cut . In Figured, the y-axis is in arbitrary units, whereas in Figure EJa the rates are 
normalised to the sum of the 3-, 4- and 5-jet rates. The renormalisation scale fi 2 has been set 
equal to q 2 , the center of mass energy of the e + e~ system. As in the previous section, the numerical 
precision is 1%. Figure 0k shows that the 5-jet rate drops drastically as y cut increases, as to be 
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Figure 4: 3-, 4- and 5-jet rates at order for the sample matrix element 



expected. The 3-jet rate decreases only slowly, as only few events are classified as 2-jet events 
and thus rejected for values of y cut < 0.08. Figure^ demonstrates how the 3-jet rate decreases 
in favour of the 4- and 5-jet rates if y cut becomes very small. 



4 Conclusions and Outlook 

In this paper, a method based on sector decomposition to calculate the double real radiation part 
of the process e + e~ — » 3 jets at O(al) has been presented. The sector decomposition algorithm 
serves to isolate, by an automated algebraic subroutine, the infrared poles which occur upon phase 
space integration if one or several particles in the final state become soft and/or collinear. In this 
way, one is dispensed from the manual construction of a subtraction scheme. The cancellation of 
the poles is shown by numerical calculation of the pole coefficients. 

For the process e + e~ — ► 3 jets at NNLO, integration over a phase space with up to five 
particles in the final state is necessary, where up to two particles can become soft and collinear. It 
has been proven that the program handles the isolation and subtraction of the poles correctly by 
considering all possible cuts of a specific diagram which is a subpart of the colour structure ~ Cp 
contained in the full matrix element at order a^. Summing over all the cuts, the poles cancel 
within the numerical precision. 

The finite part has been implemented into a Monte Carlo program which allows the inclusion 
of a measurement function in order to obtain differential distributions for arbitrary (infrared safe) 
observables. As an example, the 3-, 4- and 5-jet rates at order function of the jet 

separation parameter y cut are shown for the subpart of the full matrix element treated in this 
paper. As this toy matrix element already shows most of the problems which occur in the double 
real radiation part, while the one-loop virtual corrections combined with the 1 — > 4 phase space, 
as well as the two-loop virtual part combined with the 1 — > 3 phase space, are relatively easy 
(because the virtual integrals only lead to renormalisation factors), it is an ideal testing ground 
for the method presented here to tackle massless 1 — > 5 processes. For the calculation of the full 



matrix element, the expressions to be integrated numerically in the double real radiation part 
will of course be much larger, but the method described in this paper can handle them in the 
same way. The virtual corrections will also be more complicated, but can be treated with sector 
decomposition applied to loop integrals [27], for the kinematics of e + e _ annihilation. Therefore the 
problem is basically reduced to large file handling, book-keeping and implementation time/CPU 
time. 

For the parts of the full matrix element considered so far, the numerical stability is very good. 
A reason might be that the subtractions within the sector decomposition method are local in 
the sense of plus distributions, i.e. the singular limits in each integration variable are directly 
subtracted. 

CPU time will become an issue for the treatment of the full process, but as the method relies 
on a division of the amplitude squared into different "topologies" corresponding to different classes 
of denominator structures, the problem is naturally split into smaller subparts. If such a "trivial 
parallelisation" is not sufficient, there is still the possibility to parallelise the evaluation of the 
functions produced by sector decomposition. 

As the method is based on a universal algorithm acting on integration variables, and does not 
require analytic integration over complicated functions, it will surely see a number of interesting 
applications in the future, in particular in what concerns the production of massive particles. 
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A Massless 5-parton Phase Space 

A.l Phase space 1 — > 5 in D ^ 4 dimensions 

The phase space for the decay of one off-shell particle with momentum q into N massless particles 
with momenta p\, . . . ,pn in D dimensions is given by 



= (2^f-"f-» [ n<* i v + (pM'j-s>) 

3=1 i=l 

r N-l fifEM N-l 

= (2n) N -°^ 2 1 -" / J] d D ~ V, ^ 5+([q - £ p^ 2 ) . (A.1) 



For N = 5, we parametrise the momenta in D dimensions as (ordering of the vector components: 
(E, (D-4),x,y,z)) 

pi = £i(i,o (D - 2) ,i) 

p 2 = E 2 (l,6( D - 3 \sm6 1 ,cos6 1 ) 

P3 = E% (1,0^ D ~ 4 \ sin 02 sin 64, sin 02 cos #4, cos 62) 

Pa = E4 (1, (O ( - D ~ 5 - ) ,sin06sin05sin03),cos^6sm05sin6 , 3,cos05sin03,cos^3) 

P5 = Q ~ Pi ~ P2 - P3 - Pa (A.2) 



Inserting this parametrisation into (jA.lj) and carrying out integrations over azimuthal angles leads 
to 

J d$?_ h = (2ir f- 4D 2- 4 V(D - l)V(D - 2)V(D - 3)V(D - 4) 

4 

[sin 4 sin 5 ] D " 4 (sin 6 ) D - 5 5+ ([q - ]T K ] 2 ) (A.3) 

i=l 

where V(-D) = 2ir% /T(— ) . 
Introducing the scaled invariants j/j as new integration variables 

2/1 = s 12 /q 2 ,y2 = s 13 /q 2 ,y3 = s 23 /q 2 ,y 4 = s u /q 2 ,y 5 = s 24 /q 2 , 

y& = su/q 2 ,y7 = si 5 /q 2 ,ys = s 25 /q 2 ,y9 = W^yio = W^ 2 

leads to the Jacobian 

| det J| = 2 1( V [EiE 2 E 3 E 4 sin a sin 2 sin 3 ] 3 [sin 4 sin 5 ] 2 sin 6 . 
The Jacobian can be expressed in terms of the determinant of the Gram matrix Gij = 2pi ■ pj 

det G = —2 5 q 2 [E\E 2 E 3 E 4 sin Q\ sin 02 sin 03 sin 04 sin 05 sin 0g] 2 



=> | det J\ = yj 2 15 q 6 (- det G) \E X E 2 E 3 E 4 sin X sin 2 sin 3 ] 2 sin 4 sin 5 . 
After these variable transformations, the phase space is given by 

[d§?_+ 5 = (2TT) 5 - 4D 2- 2 - 2D V(D-l)V{D-2)V(D-3)V(D-4){q 2 ) 2D - 5 
r 10 10 

J n^i^(l-E^)(-A 5 )^" 3 e(-A 5 ) (A.4) 



j=l i=l 



where 



- A 5 = 2/102/12/22/3 + 2/I2/12/42/5 + yly2ym + 2/72/32/52/6 + 2/62/12/72/8 

+2/52/22/72/9 + 2/42/32/82/9 + 2/32/42/72/10 + 2/22/52/82/10 + 2/i2/62/92/io 
+2/10 [2/22/32/52/7 + 2/12/32/62/7 + 2/22/32/42/8 + 2/12/22/62/8 + 2/12/32/42/9 + 2/12/22/52/9] 
+2/9 [2/42/5(2/32/7 + 2/22/s) + 2/12/6(2/52/7 + 2/42/s)] + 2/62/72/8(2/32/4 + 2/22/5) 
= -idetGAg 2 ) 5 . (A.5) 

Note that V(D — 4) = 2-7r~ e /r(— e) = 0(e) is compensated by a spurious pole from (— A 5 ) _1_£ . 

For the phase space integration over the full 1 — > 5 matrix element relevant for the calculation 
of e + e~ — > 3jets at NNLO, we choose different parametrisations optimised for certain types of 
denominators occurring in the matrix element. In the following we only give the parametrisation 
which is relevant for the topology under consideration in this article. In this parametrisation, we 
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eliminate y\ by 5(1 — Ui) an d substitute ye and y 7 in favour of xe = si^/q 2 , t 7 = si^/q 2 . 

i 1 

After these substitutions, the constraint ©(— A5) is solved for 2/5, leading to 

yf = y° 5 ±VRr 5 . 

Then the condition i?5 > is solved for y$ and the condition (y% — y%) > is solved for yio- 
Making variable transformations such that all integration limits over the new variables t{ are from 
zero to one, we finally obtain 

S1345/V = h 
Sffl/? 2 = t 6 t 7 

sis/q 2 = t e t 7 (l-t 2 ) 

W^ 2 = t 3 (l-t 7 ){l-t 2 t 4 )(t e {l-t 9 )+t 9 ) 
su/q 2 = t 2 t 4 t 6 t 7 

s2A/q 2 = j/jf +ivt -yjf)*5 

S3i/q 2 = t 2 t e t 7 (1 - u) 

si5/q 2 = t 7 (1 - t 6 ) [1 - *9 (1 - tit*)] - yw 

s25/q 2 = Vs +(y£ 

= *7t9(l-t6)(l-*2*4) 

S45/V = 2/ro + {Vio ~ Via) *io (A.6) 
y± = y 8 °±d 8 /2 

y° = (l-t 6 )(l-t7){<9 + <3[t6(l-*9)-t9]}/(t6(l-*9)+t9) 

d 8 = y 8 + - % = 4 (1 - t 6 ) (1 - t 7 ) ^(1 - t 3 ) t 3 1 6 (1 - t 9 ) h/(t 6 (1 - t 9 ) + t 9 ) 

yf = y° w ±d w /2 

y% = t 2 t 7 (l-te){l-t 9 -t 4 [l-t 9 (2-t 2 )]}/(l-t 2 t 4 ) 

dio = yto ~ Vw = 4t 7 1 2 (1 - ie) ^{l-t 2 ){l-t 4 )t 4 (l-t 9 )t 9 /(l - t 2 t 4 ) . 

The solution of A5 = 0, yf, is rather lengthy and thus will not be given explicitly. 
In terms of the new variables, the phase space is given by 



[t 2 t 6 (l - * 6 )(1 - *7)] 1_2e [(l - t 2 )i 3 (l - t 3 )i 4 (l - U)t 9 (l - t 9 )]- e t 2 7 ~ 3e (A.7) 



/4 5) = (27r) 5 - 4D 2- 2 - 2Z? 2- 8e F(Z) - - 2)F(L> - 3)F(L> - 4) 



2 17 vr 9 r(-2e)r(2 - 2e) ' 
A. 2 Phase space 1 — > 5 in Z) = 4 dimensions 



In D = 4 dimensions, the Gram determinant A5 is zero due to the fact that already 4 independent 
light-like momenta pi <G {pi, • • • span Minkowski space. This leads to a nonlinear constraint 



between the Mandelstam variables in, as can be seen from eq. (|A.5[) . The momenta in D = 4 
can parametrised as in eq. (|A.2[) . but with 9q = (and no {D — 4)-dimensional component). 
The constraint 0(— A5) in eq. (|A.4|) becomes 8{— A5), leading to 2/5 = instead of y§ + (y£ — 
yjj") £5, that is, £5 takes only the values or 1. Therefore, a consistent way to obtain the 4- 
dimensional phase space from the .D-dimensional one is to integrate over t§ in l|A.7|) before any 

sector decomposition is performed, cancelling the spurious pole coming from the t$ integration 

(5) 

with V(D — 4) contained in A+ . The matrix element ME will always be of the form ME= 
Ao + Ai 2/5 + A2 y\ because in all cases where y§ is in the denominator, a different parametrisation 
will be chosen, such that the same arguments hold for a different invariant yi with i ^ 5. Using 
the fact that in our case 1 A2 = and writing y§ as 2/5 = y§ t§ + y§ (1 — £5) we obtain 

/" d$f_ 5 ME = 

r i 10 

4 5) (g 2 ) 2 ^ 5 / J] {L4 + y 5 + ^i] + [A> + y 5 "^i]} [*8(i - *a)tio(i - hoT 1 ^ 

[t 2 t e (l - i 6 )(l - t 7 )] 1_2e [(l - t 2 )t 3 (l - t 3 )U(l - t 4 )t 9 (l - t 9 )]- £ t 2 ~ 3e (A.8) 



ff (B) r(5) r(-6)r(i-6) (2^ 

r r r(l-2e) (4^)8r 2 (l/2-e)r(2-2e) 1 j 

(5) 

The new prefactor i? r ' is finite in the limit e —* 0. The matrix element only depends on the 
eight independent variables ^2, • • • ; t&, ■ ■ ■ , tio now and sector decomposition in those variables 
will isolate the "physical" infrared poles. Therefore we can still use the parametrisation HA.6j) in 
D = 4, the only difference being that S24/5 2 is given by y§ respectively y§. 

In order to construct a Monte Carlo program of (partonic) event generator type, it is useful to 
express the four momenta again in terms of angles and energies. The corresponding expressions 
in terms of Mandelstam variables are 

^ q 2 ~ S2345 ^ Q 2 ~ S1345 ^ ? 2 ~ S1245 ^ ? 2 ~ ^1235 

2q 2q ' S 2q 2q 

COS 6*1 = -1 + 2 (S1345S2345 ~ S345)/ (1 ~ S1345)/ (1 ~ «2345) 
COS 02 = -1 + 2 (Si245«2345 ~ S245) / (1 ~ S1245)/ (1 ~ S2345) 
COS6I3 = -1 + 2 (Si235«2345 - «235)/(l - Sl235)/(1 - S2345) 

The expression for cos #4 and cos #5 are more complicated and will not be given explicitly. 

Taking the limit e — > in (|A.8|) . the phase space integral over a matrix element ME in the 
above parametrisation is given by 

[t 8 (l - * 8 )*io(l " fio)]~* h te t 2 7 (1 - i 6 )(l - t 7 ) . (A.10) 



The generalisation to the case A2 7^ is trivial, leading only to additional V functions. 
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